class ASTroCAT_v16 extends DiverseEcosystem {
  
  // the verstion of the ASTroCAT model described in the 2011 paper.
  // imposed N supply, spectrum of P, spectrum of Z, and that's it. No light dependence.
  // ac_model.pde, ac_expts.pde, and ac_npz.pde from that version are combined in this file.
  
  Statevar Ptot, Ztot, N;
  Statevar[] P, Z;
  
  Flux[][] FgrazP; // one for each P-Z pair
  Flux[] Fegest; // summed over Z for each P
  Flux[] Fexcret; // summed over P for each Z
  Flux[] Fuptake, Fpmort, Fzmort;
  Flux Fsupply;
  
  FreeParam epsilon, fegest, supply; // system params
  FreeParam mu0, m, I0, zeta;  // rate params
  FreeParam ks, K0; // shape params
  FreeParam xpreyopt, dxprey; // grazing pref params
  float[][] phiP; // grazing preferences P -> Z
  boolean doAnnualCycle = false;
  
  int NP, NZ;
  
  // system definition ----------------------
  
  
  
  ASTroCAT_v16() {construct(10);}
  
  ASTroCAT_v16(int NP) {construct(NP);}
    
  void construct(int NP) {
    this.NP = NP;
    this.NZ = NP;
    
    float PMinEsd = 1; // size ranges in microns
    float PMaxEsd = 20;
    
    // state variables
    N = addVar("nutrients","N","µM");
    P = addDiverseVarSet("phytoplankton","P","µM",NP,PMinEsd, PMaxEsd);
    Z = addDiverseVarSet("zooplankton","Z","µM",NZ,1,2);
    for (int i=0; i<NP; i++) ((DiverseVar)Z[i]).esd = 2.16 * pow(((DiverseVar)P[i]).esd, 1.79);
    Ptot = addGhostVar("total P","Ptot","µM");
    Ztot = addGhostVar("total Z","Ztot","µM");
    
    // minimum biomass values: always protect a little bit of seed stock
    for (int i=0; i<NP; i++) P[i].min = 0.01;
    for (int j=0; j<NZ; j++) Z[j].min = 0.01;
    
    // fluxes
    Fsupply = addFlux("supply","supply",null,N);
    Fuptake = addFluxSet("uptake","uptake", N, P);
    FgrazP = addFluxSet2D("grazing on P","grazP",P, Z);
    Fegest = addFluxSet("egestion","egest", Z, null);
    Fexcret = addFluxSet("excretion","excret", Z, N);
    Fpmort = addFluxSet("P mortality","Pmort", P, null);
    Fzmort = addFluxSet("Z mortality","Zmort", Z, null);

    // global params
    supply = addParam("N supply rate", "supply","µM/d", 1, 0.01, 4);
    supply.logScale = true;
    supply.highlight = true;
    FreeParam p = addParam("annual cycle on/off","annCycle","", 0, 0, 1);
    p.quantize(1);
    p.highlight = true;    
    epsilon = addParam("Z growth efficiency","epsilon","", 0.33, 0, 1);
    fegest = addParam("egestion fraction","fegest","",0.33, 0, 1); // a fraction of total grazing, not of the losses
    m = addParam("P mortality as fraction of mu0","m","", 0.1, 0, 1);
    dxprey = addParam("prey size preference width","dxprey","", 0.25, 0.01, 1);
    dxprey.logScale = true;
    
    // diverse params: defined using a scale size of 1 µm for easy reading
    // std dev within a size class is defined based on the sources cited in the astrocat paper, but not used in astrocat 1.0
    //   rate params
    mu0 = addDiverseParam("P max growth","mu0","1/d",                   2.6,     0,   3,   3.8,    1,  -0.45); // initial      min      max     std dev (*,/) within a size class      allo scale size     exponent for allometry
    I0 = addDiverseParam( "Z max grazing","I0","1/d",                    26,     0,  30,   1.9,    1,  -0.4);
    zeta = addDiverseParam("Z quadratic mortality","zeta","1/µM/d",       1,     0,  10,     1,    1,   0);
    zeta.highlight = true;
    //  functional-response shape params
    ks = addDiverseParam("nutrient half-saturation","ks","µM",          0.1,   0.1,   5,     2,    1,   1);
    K0 = addDiverseParam("grazing half-saturation","K0","µM",             3,   0.3,  30,     3,    1,   0);
    K0.logScale = true;
    //  grazing preference params
    xpreyopt = addDiverseParam("optimal prey size","xpreyopt","",      0.65,   0.5,  30,     1,    1,   0.56); // has allometry but no stochasticity; noise in prey prefs is handled separately

    phiP = new float[NP][NZ]; // prey preferences: internal array
    addParam("max preference","maxPref","", 10, 0, 10);

    p = addParam("prefs allometry on/off","prefsAllo","",1, 0, 1);
    p.quantize(1);
    p.highlight = true;
    p = addParam("other allometry on/off","otherAllo","",1, 0, 1);
    p.quantize(1);
    p.highlight = true;    
    p = addParam("prefs stochasticity","prefsRand","",0, 0, 1); // other stochasticity (in mu0, I0, for example) is omitted in this version of the model
    p.quantize(1);
    p.highlight = true;        
 
   // set a better scale size, so the monotype version doesn't run with 1 µm cells
    float xP = sqrt(PMinEsd*PMaxEsd); // geometric means of the size ranges
    float xZ = 2.16 * pow(xP, 1.79); // set Z sizes such that each P is optimal prey for one Z
    ((DiverseParam)mu0).initial *= pow(xP / ((DiverseParam)mu0).alloScaleSize, ((DiverseParam)mu0).alloExp);
    ((DiverseParam)mu0).alloScaleSize = xP;
    ((DiverseParam)I0).initial *= pow(xZ / ((DiverseParam)I0).alloScaleSize,  ((DiverseParam)I0).alloExp);
    ((DiverseParam)I0).alloScaleSize = xZ;
    ((DiverseParam)xpreyopt).initial *= pow(xZ / ((DiverseParam)xpreyopt).alloScaleSize,  ((DiverseParam)xpreyopt).alloExp);
    ((DiverseParam)xpreyopt).alloScaleSize = xZ;
    ((DiverseParam)ks).initial *= pow(xP / ((DiverseParam)ks).alloScaleSize,  ((DiverseParam)ks).alloExp);
    ((DiverseParam)ks).alloScaleSize = xP;
    // adjust the names of the allometric params so the user can see what the sliders control
    mu0.name += " at " + round(xP) + " µM";
    I0.name += " at " + round(xZ) + " µM";
    ks.name += " at " + round(xP) + " µM";
    xpreyopt.name += " at " + round(xZ) + " µM";
  }
  
  
  
  void calcFluxes() {
    for (int k=0; k<Ncells; k++) { // for each spatial cell...
      // supply
      Fsupply.current[k] = supply.current;
      if (doAnnualCycle) Fsupply.current[k] *= 0.5*(1+cos(time.current[0] * TWO_PI/365.));
      // phyto growth
      for (int i=0; i<NP; i++) {
        float mu0i = getCurrent(mu0, P[i]);
        float ksi = getCurrent(ks, P[i]);
        float nutlim = N.current[k] / (ksi + N.current[k]);
        Fuptake[i].current[k] = mu0i * nutlim * P[i].current[k];
      }
      
      // grazing
      //     assemble all diet elements
      float[][] PscaledAsFood = new float[NP][NZ];
      float[] totalFood = new float[NZ];
      for (int j=0; j<NZ; j++) {
        float K0j = getCurrent(K0, Z[j]);
        totalFood[j] = 0;
        for (int i=0; i<NP; i++) {
          PscaledAsFood[i][j] = phiP[i][j] / K0j * P[i].current[k];
          totalFood[j] += PscaledAsFood[i][j];
        }
      }
      //     initialize excretion and egestion fields
      for (int i=0; i<NP; i++) Fegest[i].current[k] = 0;
      for (int j=0; j<NP; j++) Fexcret[j].current[k] = 0;
      //     calculate fluxes
      for (int j=0; j<NZ; j++) {
        float I0j = getCurrent(I0, Z[j]);
        for (int i=0; i<NP; i++) {
          FgrazP[i][j].current[k] = I0j * Z[j].current[k] * PscaledAsFood[i][j] / (1. + totalFood[j]);
          Fegest[j].current[k] += fegest.current * FgrazP[i][j].current[k];
          Fexcret[j].current[k] += (1.0 - fegest.current - epsilon.current) * FgrazP[i][j].current[k];          
        }
      }
      
      // mortality
      for (int i=0; i<NP; i++) {
        float mi = getCurrent(m, P[i]);
        float mu0i = getCurrent(mu0, P[i]);
        Fpmort[i].current[k] = mi * mu0i * P[i].current[k];
//        Fpmort[i].current[k] += 0.03 * supply.current * P[i].current[k]; // chemostat loss term
      }
      for (int j=0; j<NZ; j++) {
        float zetaj = getCurrent(zeta, Z[j]);
        Fzmort[j].current[k] = zetaj * Z[j].current[k] * Ztot.current[k];
//        Fzmort[j].current[k] += 0.03 * supply.current * Z[j].current[k]; // chemostat loss term
      }
      
      // diagnostics
      Ptot.current[k] = 0;
      for (int i=0; i<NP; i++) Ptot.current[k] += P[i].current[k];
      Ztot.current[k] = 0;
      for (int j=0; j<NZ; j++) Ztot.current[k] += Z[j].current[k];
    }
    
  }



  void updatePrivateParams() {
    doAnnualCycle = (getParamVal("annCycle")==1);
    
    // turn off stochasticity in rates and half-saturations
    ((DiverseParam)mu0).randMultiplier = 0;
    ((DiverseParam)ks).randMultiplier = 0;
    ((DiverseParam)I0).randMultiplier = 0;
    ((DiverseParam)K0).randMultiplier = 0;
    ((DiverseParam)zeta).randMultiplier = 0;
    
    // turn allometry on and off, in response to the user's choice
    float r = getParamVal("otherAllo");
    ((DiverseParam)mu0).doAllo = (r==1);
    ((DiverseParam)ks).doAllo = (r==1);
    ((DiverseParam)I0).doAllo = (r==1);
    ((DiverseParam)K0).doAllo = (r==1);
    ((DiverseParam)zeta).doAllo = (r==1);
    
    // set grazing preferences, including allometry and stochasticity
    r = getParamVal("prefsRand");
    boolean doAllo = (getParamVal("prefsAllo")==1);
    float maxPref = getParam("maxPref").current;
    for (int j=0; j<NZ; j++) {
      float xpreyoptj = getCurrent(xpreyopt, Z[j]);
      float dxpreyj = getCurrent(dxprey, Z[j]);
      for (int i=0; i<NP; i++) {
        float xpreyi = ((DiverseVar)P[i]).esd;
        float R = randn(((DiverseVar)P[i]).randSignature,((DiverseVar)Z[j]).randSignature);
        phiP[i][j] = exp(r*R);
        if (doAllo) phiP[i][j] *= exp(-sq((log10(xpreyi)-log10(xpreyoptj)) / dxpreyj));
        phiP[i][j] = constrain(phiP[i][j], 0, maxPref);
      }
    }
    
  }

}



class ASTroCAT_v16_NetworkView extends DiverseNetworkView {
  
  ASTroCAT_v16 eco;
  
  ASTroCAT_v16_NetworkView(Ecosystem eco) {
    this.eco = (ASTroCAT_v16)eco;
  }
  
  void define() {
    varScale = 10 * (40.0/eco.NP);
    fluxScale = 0.5 * (40.0/eco.NP);

    float R = 0.8;
    float xcen = 0.2;
    float ycen = 0.1;
    int margin = 5;

    float[] px = new float[eco.NP];
    float[] py = new float[eco.NP];
    float[] zx = new float[eco.NZ];
    float[] zy = new float[eco.NZ];
      
    float[][] vertices = fluxArc(xcen-R, ycen-R, xcen-R+2*R*cos(PI/8), ycen-R+2*R*sin(PI/8), 30, eco.NZ+margin);
    for (int j=0; j<eco.NZ; j++) {
      zx[j] = vertices[j+margin][0];
      zy[j] = vertices[j+margin][1];
      Dot d = addDot(eco.Z[j], zx[j], zy[j], transparency(styles.dkBrownColor,0.8));
      d.name = nicenum(round(((DiverseVar)eco.Z[j]).esd));
    }
    vertices = fluxArc(xcen-R, ycen-R, xcen-R+2*R*cos(HALF_PI-PI/8), ycen-R+2*R*sin(HALF_PI-PI/8), 30, eco.NP+margin);
    for (int i=0; i<eco.NP; i++) {
      px[i] = vertices[i+margin][0];
      py[i] = vertices[i+margin][1];
      Dot d = addDot(eco.P[i], px[i], py[i], transparency(styles.medGreenColor,0.8));
      d.name = nicenum(round(((DiverseVar)eco.P[i]).esd));
    }
    addDot(eco.N, xcen-R, ycen-R, styles.dkBlueColor);
    
    // supply arrow
    Arrow a = addArrow(eco.Fsupply, xcen-2*R, ycen-R, styles.dkBlueColor, 0);
    a.arrowheadPos = 0.8;

    // regeneration and loss arrows
    addArrowSet(eco.Fexcret, styles.ltGrayColor, 150);
    for (int j=0; j<eco.NZ; j++) {
      a = addArrow(eco.Fzmort[j], zx[j]+0.4, zy[j], styles.medGrayColor, 0);
      a.setFrom(zx[j], zy[j]);
    }
    for (int i=0; i<eco.NP; i++) {
      a = addArrow(eco.Fpmort[i], px[i]+0.4, py[i], styles.medGrayColor, 0);
      a.setFrom(px[i], py[i]);
    }
    
    // growth & grazing arrows
    Arrow[] aa = addArrowSet(eco.Fuptake, styles.ltGreenColor, 120);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    aa = addArrowSet(eco.FgrazP, styles.dkBrownColor, 90);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
  }
  
  void updateArrows() {
    for (int i=0; i<arrows.length; i++) {
      arrows[i].update(fluxes[i].current[currentCell]);
    }
    if (arrows.length > 400) {
      arrows[0].draw(); // draw the supply arrow only if there's a lot of arrows
    } else {
      for (int i=0; i<arrows.length; i++) arrows[i].draw();
    }
  }
  
}



// --------------------------------------------------------------------------------------------------------------------------




// ------------------------------------------------------------------------------
// chemostat environment to run ASTroCAT experiments in -------------------------

class ASTroCAT_v16_Chemostat extends Environment {
  float seed = 0.01;
  float Ntot = 0;
  
  ASTroCAT_v16_Chemostat(int NP) {
    linkToEcosystem(new ASTroCAT_v16(NP));
    for (int j=0; j<NP; j++) {
      getVar("P-" + j).setInitial(seed/NP);
      getVar("Z-" + j).setInitial(seed/NP);
    }
    getVar("N").setInitial(max(Ntot-2*seed, 0));
    startArchive(1);
    archive.saveFluxes = false;
  }

  void saveArchive() {
    if (archive != null) {
      ASTroCAT_v16 eco = (ASTroCAT_v16)this.eco;
      // save rates and sizes by size class
      archive.addNameValuePair("params-mu0", eco.getCurrent(eco.mu0, eco.P));
      archive.addNameValuePair("params-I0", eco.getCurrent(eco.I0, eco.Z));
      archive.addNameValuePair("params-zeta", eco.getCurrent(eco.zeta, eco.Z));
      archive.addNameValuePair("params-ks", eco.getCurrent(eco.ks, eco.P));
      archive.addNameValuePair("params-K0", eco.getCurrent(eco.K0, eco.Z));      
      archive.addNameValuePair("phiP", eco.phiP); 
      float[] Pesd = new float[eco.NP];
      for (int i=0; i<Pesd.length; i++) Pesd[i] = ((DiverseVar)eco.P[i]).esd;
      float[] Zesd = new float[eco.NP];
      for (int i=0; i<Zesd.length; i++) Zesd[i] = ((DiverseVar)eco.Z[i]).esd;
      archive.addNameValuePair("Pesd", Pesd);
      archive.addNameValuePair("Zesd", Zesd);
      // calculate time-mean mu, g for each P class. This is so that it isn't necessary to save the
      // fluxes to the netcdf file (the 1600 grazing fluxes take too long to save!)
      int N = archive.states.length;
      float[] mu = new float[eco.NP];
      float[] g = new float[eco.NP];
      for (int n=0; n<N; n++) {
        for (int i=0; i<eco.NP; i++) {
          mu[i] += eco.Fuptake[i].current[0] / eco.P[i].current[0];
          for (int j=0; j<eco.NP; j++) {
            g[i] += eco.FgrazP[i][j].current[0] / eco.P[i].current[0];
          }
        }
      }
      for (int i=0; i<eco.NP; i++) {
        mu[i] /= N;
        g[i] /= N;
      }
      archive.addNameValuePair("diag-mu", mu);
      archive.addNameValuePair("diag-g", g);
    }
    super.saveArchive();
  }
  
}




// ------------------------------------------------------------------------------
// scripted experiments ---------------------------------------------------------


class ASTroCAT_Case extends Script {
  float intLength = 3650;
  float saveTimestep = 1;
  int NP = 40;
  String name;
  Environment env;
  
  ASTroCAT_Case(String casename) {
    name = casename;
    env = new ASTroCAT_v16_Chemostat(NP);
    initialize();
  }
  
  void initialize() {
    env.resetParams();
    env.getParam("otherAllo").current = 1;
    env.getParam("prefsAllo").current = 1;        
    env.getParam("prefsRand").current = 0;
    env.getParam("maxPref").current = 10;
  }
  
  void reset() {
    env.updatePrivateParams();
    env.rewind();
    if (env.archive==null) env.startArchive(saveTimestep);
    env.archive.saveFluxes = false;
    env.archive.basename = name;
  }
  
  void run() {
    reset();
    println("running " + name);
    while (env.time.current[0] < intLength) {
      env.integrate(env.archive.dt);
    }
    env.saveArchive();
  }  
}


class ASTroCAT_InterannualCase extends ASTroCAT_Case {
  int nYears = 100;
  float supply_mean = 1;
  float supply_std = 0.25;
  
  ASTroCAT_InterannualCase(String casename) {
    super(casename);
  }
  
  void initialize() {
    super.initialize();
    env.getParam("annCycle").current = 1;
    env.getParam("supply").current = supply_mean;
  }
  
  void run() {
    println("running " + name);
    reset();
    intLength = 365./2.;
    while (env.time.current[0] < intLength) env.integrate(env.archive.dt);
    for (int i=0; i<nYears; i++) {
      intLength = 365 * (i+0.5);
      env.getParam("supply").current = supply_mean + supply_std * randn();
      while (env.time.current[0] < intLength) env.integrate(env.archive.dt);
      println("year " + i);
      if (i % 10 == 0) {println("saving"); env.saveArchive();}
    }
    env.saveArchive();
  }
}




void ASTroCAT_v16_allExpts() {
  ASTroCAT_expt1_steady();
  ASTroCAT_expt2_steady_varying_dxprey();
  ASTroCAT_expt3_steady_varying_supply();
  ASTroCAT_expt4_steady_varying_zeta();
  ASTroCAT_expt5_steady_stochastic();
  ASTroCAT_expt6_anncycle_varying_dxprey();
  ASTroCAT_expt7_interannual_varying_dxprey();
  ASTroCAT_expt8_interannual_varying_NP();
  ASTroCAT_expt9_interannual_stochastic();
}




void ASTroCAT_expt1_steady() {
  ASTroCAT_Case S = new ASTroCAT_Case("steady-basecase");
  S.run();
  S.name = "steady-prefsAlloOnly";
  S.env.getParam("otherAllo").current = 0;
  S.run();
}

void ASTroCAT_expt2_steady_varying_dxprey() {
  float[] widths = {0.063, 0.088, 0.125, 0.177, 0.25, 0.354, 0.5};
  ASTroCAT_Case S;
  for (int i=0; i<widths.length; i++) {
    S = new ASTroCAT_Case("steady-dxprey-" + i);
    S.env.getParam("dxprey").current = widths[i];
    S.run();
  }  
}

void ASTroCAT_expt3_steady_varying_supply() {
  float[] levels = {0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4};
  ASTroCAT_Case S;
  for (int i=0; i<levels.length; i++) {
    S = new ASTroCAT_Case("steady-supply-" + i);
    S.env.getParam("supply").current = levels[i];
    S.run();
  }  
}

void ASTroCAT_expt4_steady_varying_zeta() {
  float[] zetas = {0.25, 0.5, 1, 2, 4};
  ASTroCAT_Case S;
  for (int i=0; i<zetas.length; i++) {
    S = new ASTroCAT_Case("steady-zeta-" + i);
    S.env.getParam("zeta").current = zetas[i];
    S.run();
  }  
}

void ASTroCAT_expt5_steady_stochastic() {
  float[] levels = {0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4};
  int Nreps = 8;
  ASTroCAT_Case S = new ASTroCAT_Case("");
  for (int j=0; j<Nreps; j++) { // replicates
    randomSeed(round(1e4*j));
    ((DiverseEcosystem)S.env.eco).reseed();
    S.env.getParam("supply").current = 1;
    S.name = "steady-stoch-" + j + "-r-01";
    S.env.getParam("prefsRand").current = 0.1;
    S.run();    
    S.name = "steady-stoch-" + j + "-r-025";
    S.env.getParam("prefsRand").current = 0.25;
    S.run();
    for (int i=0; i<levels.length; i++) {
      S.name = "steady-stoch-" + j + "r-1-supply-" + i;
      S.env.getParam("prefsRand").current = 1;
      S.env.getParam("supply").current = levels[i];
      S.run();
    }
  }
} 

void ASTroCAT_expt6_anncycle_varying_dxprey() {
  float[] widths = {0.063, 0.088, 0.125, 0.177, 0.25, 0.354, 0.5};
  ASTroCAT_Case S;
  for (int i=0; i<widths.length; i++) {
    S = new ASTroCAT_Case("anncyc-dxprey-" + i);
    S.env.getParam("annCycle").current = 1;
    S.env.getParam("dxprey").current = widths[i];
    S.intLength = 365*100.5;
    S.run();
  }  
}

void ASTroCAT_expt7_interannual_varying_dxprey() {
  float[] widths = {0.063, 0.088, 0.125, 0.177, 0.25, 0.354, 0.421, 0.5, 0.594};
  ASTroCAT_InterannualCase S;
  for (int i=0; i<widths.length; i++) {
  S = new ASTroCAT_InterannualCase("interann-dxprey-" + i);
    S.env.getParam("dxprey").current = widths[i];
    S.run();
  }  
}

void ASTroCAT_expt8_interannual_varying_NP() {
//  int[] NPs = {2, 4, 8, 16, 24, 40, 64};
  int[] NPs = {3, 6};
  ASTroCAT_InterannualCase S;
  for (int i=0; i<NPs.length; i++) {
    S = new ASTroCAT_InterannualCase("interann-NP-" + NPs[i]);
    S.env = new ASTroCAT_v16_Chemostat(NPs[i]);
    S.initialize();
    S.run();
  }  
}

void ASTroCAT_expt9_interannual_stochastic() {
  ASTroCAT_InterannualCase S = new ASTroCAT_InterannualCase("interann-stoch-01");
  S.env.getParam("prefsRand").current = 0.1;
  S.run();
  S.name = "interann-stoch-025";
  S.env.getParam("prefsRand").current = 0.25;
  S.run();
  S.name = "interann-stoch-1";
  S.env.getParam("prefsRand").current = 1;
  S.run();
}



// -------------------------------------------------------------------------------------------------------------------------



// simple, 1-box NPZ and NPZD models, no light dependence, no sinking, no physics.
// included as examples with ASTroCAT.

class SimpleNPZ extends Ecosystem {

  Statevar N,P,Z,totalN; 
  Flux Fuptake, Fingest, Fpmort, Fzmort, Fexcret;
  FreeParam mu0, I0, ks, Ks, m, mz, epsilon;
  
  // system definition ----------------------
  
  SimpleNPZ() {
    N = addVar("nutrients","N","µM");
    P = addVar("phytoplankton","P","µM");
    P.sinks = true;
    Z = addVar("zooplankton","Z","µM");
    totalN = addGhostVar("total nitrogen","T","µM");
    
    Fuptake = addFlux("uptake","uptake", N, P);
    Fingest = addFlux("ingestion","ingest", P, Z);
    Fpmort = addFlux("P mortality","Pmort", P, N);
    Fzmort = addFlux("Z mortality","Zmort", Z, N);
    Fexcret = addFlux("excretion","excret", Z, N);

    mu0 = addParam("P max growth","mu0","1/d",2,0,2);
    ks = addParam("nutrient half-saturation","ks","µM",1, 0.1, 5);
    m = addParam("P mortality","m","1/d",0.1, 0, 0.4);
    I0 = addParam("Z max growth","I0","1/d", 3, 0, 5);
    Ks = addParam("grazing half-saturation","Ks","µM",3, 0, 4);
    epsilon = addParam("Z growth efficiency","epsilon","",0.3, 0, 1);
    mz = addParam("Z mortality","mz","1/d",0.4, 0, 1.0);
  }
  
  // grazing and mortality are linear: simple, prone to oscillations
  void calcFluxes() {
    for (int i=0; i<Ncells; i++) {
      Fuptake.current[i] = mu0.current * N.current[i] / (ks.current + N.current[i]) * P.current[i];
      Fingest.current[i] = I0.current * P.current[i] / (Ks.current + P.current[i]) * Z.current[i];
      Fpmort.current[i] = m.current * P.current[i];
      Fzmort.current[i] = mz.current * Z.current[i];
      Fexcret.current[i] = (1 - epsilon.current) * Fingest.current[i];
      totalN.current[i] = N.current[i] + P.current[i] + Z.current[i];
    }
  }
  
}


class SimpleNPZ_NetworkView extends NetworkView {
  
  SimpleNPZ eco;
  
  SimpleNPZ_NetworkView(Ecosystem eco) {
    this.eco = (SimpleNPZ)eco;
  }

  void define() {
    varScale = 50;
    fluxScale = 4;
    float theta0 = -120;
    addDot_polar(eco.Z, 0.7, theta0-240, styles.medRedColor);
    addDot_polar(eco.P, 0.7, theta0-120, styles.ltGreenColor);
    addDot_polar(eco.N, 0.7, theta0, styles.dkBlueColor);
    
    addArrow(eco.Fuptake, styles.ltGreenColor, 120);
    addArrow(eco.Fingest, styles.medRedColor, 120);
    addArrow(eco.Fexcret, styles.medBlueColor, 120);
    addArrow(eco.Fpmort, styles.medBlueColor, 120);
    addArrow(eco.Fzmort, styles.medBlueColor, 120);
  }
      
}




// ---------------------------------------------------------------------------------------------------------




class SimpleNPZD extends Ecosystem {
  
  Statevar N,P,Z,D,totalN; 
  Flux Fuptake, Fingest, Fegest, Fpmort, Fzmort, Fexcret, Fremin, Fsupply, Fsink;
  FreeParam S, mu0, ks, alpha, attSW, attP, m, I0, Ks, epsilon, fegest, zeta, Rremin, wsinkD;
  float nominalBoxDepth = 20; // used in converting the D sinking velocity to a fraction lost per day
  
  // system definition ----------------------
  
  SimpleNPZD() {
    N = addVar("nutrients","N","µM");
    P = addVar("phytoplankton","P","µM");
    Z = addVar("zooplankton","Z","µM");
    D = addVar("detritus","D","µM");
    D.sinks = true;
    totalN = addGhostVar("total nitrogen","T","µM");

    Fsupply = addFlux("N supply","supply", null, N);
    Fuptake = addFlux("uptake","uptake", N, P);
    Fingest = addFlux("ingestion","ingest", P, Z);
    Fegest = addFlux("egestion","egest", P, D);
    Fexcret = addFlux("excretion","excret", Z, N);
    Fpmort = addFlux("P mortality","Pmort", P, D);
    Fzmort = addFlux("Z mortality","Zmort", Z, D);
    Fremin = addFlux("remineralization","remin", D, N);
    Fsink = addFlux("D sinking","Dsink",D,null);

    S = addParam("N supply rate","S","µM/d",1,0,2);
    S.highlight = true;
    mu0 = addParam("P max growth","mu0","1/d",2,0,3);
    ks = addParam("nutrient half-saturation","ks","µM",1, 0.1, 5);
    m = addParam("P mortality","m","1/d",0.1, 0, 0.4);
    I0 = addParam("Z max growth","I0","1/d", 3, 0, 6);
    Ks = addParam("grazing half-saturation","Ks","µM",3, 0.1, 10);
    Ks.logScale = true;
    epsilon = addParam("Z growth efficiency","epsilon","", 0.3, 0, 1);
    fegest = addParam("egestion fraction","fegest","",0.5, 0, 1);
    zeta = addParam("Z quadratic mortality","zeta","1/µM/d",2.0, 0, 4.0);
    Rremin = addParam("remineralization rate","remin","1/d",0.1, 0.01, 1);
    Rremin.logScale = true;
    wsinkD = addParam("D sinking rate","wsinkD","m/day",10,0,100);
  }
  
  void calcFluxes() {
    for (int i=0; i<Ncells; i++) {
      Fsupply.current[i] = S.current;
      float nutlim = N.current[i] / (ks.current + N.current[i]);
      Fuptake.current[i] = mu0.current * nutlim * P.current[i];
      float Fgraz = I0.current * sq(P.current[i]) / (sq(Ks.current) + sq(P.current[i])) * Z.current[i];
      Fingest.current[i] = epsilon.current * Fgraz;
      Fegest.current[i] = (1 - epsilon.current) * fegest.current * Fgraz;
      Fexcret.current[i] = (1 - epsilon.current) * (1 - fegest.current) * Fgraz;
      Fingest.current[i] += Fexcret.current[i]; // because excretion is written as a flux from Z to N
      Fpmort.current[i] = m.current * P.current[i];
      Fzmort.current[i] = zeta.current * sq(Z.current[i]);
      Fremin.current[i] = Rremin.current * D.current[i];
      Fsink.current[i] = wsinkD.current / nominalBoxDepth * D.current[i]; // pretty rough representation of sinking
      totalN.current[i] = N.current[i] + P.current[i] + Z.current[i] + D.current[i];      
    }
  }

}

class SimpleNPZD_NetworkView extends NetworkView {
  
  SimpleNPZD eco;
  
  SimpleNPZD_NetworkView(Ecosystem eco) {
    this.eco = (SimpleNPZD)eco;
  }

  void define() {
    varScale = 50;
    fluxScale = 4;
    float theta0 = -120;
    Dot Ddot = addDot_polar(eco.D, 0.7, theta0-270, styles.ltBlueColor);
    addDot_polar(eco.Z, 0.7, theta0-180, styles.medRedColor);
    addDot_polar(eco.P, 0.7, theta0-90, styles.ltGreenColor);
    Dot Ndot = addDot_polar(eco.N, 0.7, theta0, styles.dkBlueColor);
    
    addArrow(eco.Fsupply, Ndot.xcen, Ndot.ycen-0.3, styles.dkBlueColor, 0);
    addArrow(eco.Fuptake, styles.ltGreenColor, 90);
    addArrow(eco.Fingest, styles.medRedColor, 90);
    addArrow(eco.Fegest, styles.medBlueColor, 30);
    addArrow(eco.Fexcret, styles.medBlueColor, 30);
    addArrow(eco.Fpmort, styles.medBlueColor, -30);
    addArrow(eco.Fzmort, styles.medBlueColor, 90);
    addArrow(eco.Fremin, styles.dkBlueColor, 90);
    addArrow(eco.Fsink, Ddot.xcen, Ddot.ycen-0.3, styles.medBlueColor, 0);
  }
  
}



// ----------------------------------------------------------------------------------------------------------------------------------


class DiverseEcosystem extends Ecosystem {
  // utilities for making an allometrically or stochastically diverse ecosystem.
  // ASTroCAT 1.6 doesn't take advantage of the stochasticity part (except for prey prefs, where stochasticity is introduced by another mechanism anyway).
  // probably this generalization layer isn't worth it--each version of ASTroCAT should just incorporate what it needs from this in itself.
  
  
  DiverseEcosystem() {}

  // arrays of diverse statevars
  
  DiverseVar[] addDiverseVarSet(String name, String shortname, String units, int N, float minEsd, float maxEsd) {
    DiverseVar[] V = new DiverseVar[N];
    for (int i=0; i<N; i++) {
      float esd = exp(map(i, 0, N-1, log(minEsd), log(maxEsd)));
      V[i] = new DiverseVar(name + "-" + i, shortname + "-" + i, units, Ncells, esd);
      V[i].advects = true;
      V[i].mixes = true;
      V[i].sinks = false;
      V[i].min = 0;
      vars = (Statevar[])append(vars,V[i]);
      V[i].hidden = false;
    }
    return V;
  }

  // arrays of fluxes

  Flux[] addFluxSet(String name, String shortname, Statevar[] from, Statevar to) {
    Flux[] F = new Flux[from.length];
    for (int i=0; i<F.length; i++) F[i] = addFlux(name + "-" + i, shortname + "-" + i, from[i], to);
    return F;
  }

  Flux[] addFluxSet(String name, String shortname, Statevar from, Statevar[] to) {
    Flux[] F = new Flux[to.length];
    for (int i=0; i<F.length; i++) F[i] = addFlux(name + "-" + i, shortname + "-" + i, from, to[i]);
    return F;    
  }

  Flux[][] addFluxSet2D(String name, String shortname, Statevar from[], Statevar[] to) {
    Flux[][] F = new Flux[from.length][to.length];
    for (int i=0; i<from.length; i++) {
      for (int j=0; j<to.length; j++) {
        F[i][j] = addFlux(name + "-" + i + "-" + j, shortname + "-" + i + "-" + j, from[i], to[j]);
      }
    }
    return F;    
  }

  // diverse params
  
  DiverseParam addDiverseParam(String name, String shortname, String units, float initial, float min, float max, float randLogStd) {
    return addDiverseParam(name,shortname,units,initial,min,max,randLogStd,1,0);
  }
  
  DiverseParam addDiverseParam(String name, String shortname, String units, float initial, float min, float max, float randLogStd, float alloScaleSize, float alloExp) {
    DiverseParam P = new DiverseParam(name,shortname,units,initial,min,max,randLogStd,alloScaleSize,alloExp);
    params = (FreeParam[])append(params,P);
    return P;
  }

  float getCurrent(FreeParam P, Statevar V) {
    if (P instanceof DiverseParam && V instanceof DiverseVar) {
      return getCurrent((DiverseParam)P, (DiverseVar)V);
    } else {
      return P.current;
    }
  }
  
  float getCurrent(DiverseParam P, DiverseVar V) { // gets the current value of P appropriate to the statevar V, taking into account random variation and allometry
    float val = P.current; // no allometry, no randomness yet
    if (P.doAllo) val *= pow(V.esd / P.alloScaleSize, P.alloExp); // apply allometry
    val *= pow(P.randLogStd, P.randMultiplier * randn(P.randSignature,V.randSignature)); // apply randomness
    val = constrain(val, P.min, P.max);
    return val;
  }
  
  float[] getCurrent(FreeParam P, Statevar[] V) {
    float[] val = new float[V.length];
    for (int i=0; i<val.length; i++) val[i] = getCurrent(P, V[i]);
    return val;
  }
  
  void reseed() {
    for (int i=0; i<vars.length; i++) {
      if (vars[i] instanceof DiverseVar) ((DiverseVar)vars[i]).reseed();
    }
    for (int i=0; i<params.length; i++) {
      if (params[i] instanceof DiverseParam) ((DiverseParam)params[i]).reseed();
    }
  }
      
}

// ------------------------------------------------

class DiverseVar extends Statevar { // stores a size (esd) and a random signature
  
  float esd = 1;
  float randSignature = random(1);
  
  DiverseVar(String name, String shortname, String units, float esd) {
    this.name = name;
    this.shortname = shortname;
    this.units = units;
    this.esd = esd;
  }

  DiverseVar(String name, String shortname, String units, int N, float esd) {
    this.name = name;
    this.shortname = shortname;
    this.units = units;
    allocate(N);
    this.esd = esd;
  }
  
  void reseed() {
    randSignature = random(1);
  }
    
}

// ------------------------------------------------

class DiverseParam extends FreeParam { // stores a random signature and a specification of an allometric curve
  
  float randSignature = random(1);
  float randLogStd = 0; // std dev of the random distribution of log values around the allometric central tendency, in data units
  // min and max are constraints on allowable values (e.g., min=0 to keep growth rates from going negative), not related to the scale for random variation
  float alloScaleSize = 1; // for allometry, min max initial current are all values at this ESD
  float alloExp = 0; // exponent in the allometric power law:    rate = (ESD/alloScaleSize)^alloExp
  
  float randMultiplier = 1; // multiplies randLogStd. Use this to turn stochasticity down or off without losing the record of the intrinsic scale of variation
  boolean doAllo = true;
  
  DiverseParam(String name, String shortname, String units, float initial, float min, float max, float randLogStd) {
    construct(name, shortname, units, initial, min, max);
    this.randLogStd = randLogStd;
  }

  DiverseParam(String name, String shortname, String units, float initial, float min, float max, float randLogStd, float alloScaleSize, float alloExp) {
    construct(name, shortname, units, initial, min, max);
    this.randLogStd = randLogStd;
    this.alloScaleSize = alloScaleSize;
    this.alloExp = alloExp;
  }
  
  void reseed() {
    randSignature = random(1);
  }
  
}

// ------------------------------------------------

class DiverseNetworkView extends NetworkView {
  
  DiverseNetworkView() {}  

  Dot[] addDotSet_polar(Statevar[] V, float Rmin, float Rmax, float theta, color col) {
    Dot[] D = new Dot[V.length];
    for (int i=0; i<V.length; i++) {
      D[i] = addDot_polar(V[i], map(i, 0, V.length-1, Rmin, Rmax), theta, col);
    }
    return D;
  }

  Arrow[] addArrowSet(Flux[] F, color col, float angle) {
    Arrow[] A = new Arrow[F.length];
    for (int i=0; i<F.length; i++) {
      A[i] = addArrow(F[i], col, angle);
    }
    return A;
  }
  
  Arrow[] addArrowSet(Flux[] F, float dx, float dy, color col, float angle) {
    Arrow[] A = new Arrow[F.length];
    for (int i=0; i<F.length; i++) {
      A[i] = addArrow(F[i], dx, dy, col, angle);
    }
    return A;
  }

  Arrow[] addArrowSet(Flux[][] F, color col, float angle) {
    Arrow[] A = new Arrow[F.length * F[0].length];
    int ind = 0;
    for (int i=0; i<F.length; i++) {
      for (int j=0; j<F[0].length; j++) {
        A[ind] = addArrow(F[i][j], col, angle);
        ind++;
      }
    }
    return A;
  }
  
}

